Nonequilibrium work distribution of a quantum harmonic oscillator 
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We analytically calculate the work distribution of a quantum harmonic oscillator with arbitrary 
time-dependent angular frequency. We provide detailed expressions for the work probability density 
for adiabatic and nonadiabatic processes, in the limit of low and high temperature. We further verify 
the validity of the quantum Jarzynski equality. 
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I. INTRODUCTION 

Fluctuations play a fundamental role in small nonequi- 
librium systems. While thermodynamic quantities like 
the free energy difference can easily be calculated for 
slow, quasistatic processes from the knowledge of the 
initial and final equilibrium partition functions, A_F — 
Fi — Fq = —kTln Zi/Zq, a general method for evaluat- 
ing AF in the case of fast, nonequilibrium transforma- 
tions was lacking until recently. A decade ago, Jarzynski 
introduced a remarkable equality that relates the free en- 
ergy difference to the averaged exponentiated work done 
during the process 



AF = -fcrin(e 



(1) 



with (e"'^^) = / dWe-'^^ P(W) and (3 = 1/kT the in- 
verse temperature. Equation ^ allows to compute the 
difference AF for arbitrary transformations, quasistatic 
or not, once the corresponding distribution of work P(W) 
is known. The free energy difference is thus entirely 
determined by the fluctuations of the total work W, a 
property which emphasizes the importance of the work 
probability density. The usefulness of Eq. ([1]) for evalu- 
ating A_F in single biomolecules has been demonstrated 
in DNA unfolding experiments [2, 3]. 

Extensions of the classical Jarzynski equality ^ to 
quantum-mechanical systems have been discussed in 
Refs. 4, 5, 6, 7, 8, 9, 10]. In that respect, it is worth 
mentioning that the concept of work in quantum ther- 
modynamics has been clarified in Ref. |9|. Thus for a 
quantum Hamiltonian system, work cannot be consid- 
ered as the expectation value of some work operator, in 
contrast to some other thermodynamic quantities like en- 
ergy. Rather the exponential average of the work, as well 
as its characteristic function, are given by time-ordered 
quantum correlation functions. 

Our aim in this paper is to determine the work dis- 
tribution for a simple quantum system, namely a time- 
dependent harmonic oscillator. We obtain the exact an- 
alytical expression of the characteristic function of the 
work and provide approximations for the work distribu- 
tion for adiabatic and nonadiabatic processes, in the limit 
of low and high temperature. Moreover, we check the va- 
lidity of the Jarzynski equality for arbitrary processes. 



II. WORK CHARACTERISTIC FUNCTION 

We consider a quantum-mechanical harmonic oscilla- 
tor with a time-dependent angular frequency The 
corresponding Hamiltonian is of the usual form 



2m 2 



(2) 



We shall be interested in calculating the work distri- 
bution of the oscillator when the angular frequency is 
changed from an initial value at t = to a final value 
uoi a,t t = T. We denote by (j)*^ the eigenf unctions and by 
E*^ = ?iw(t)(n-t-i) the eigenvalues of the Hamiltonian ^ 
at any given time t. We make the additional assumption 
that the oscillator is initially thermalized with inverse 
temperature (3 but is otherwise isolated. 

The probability density of the total work done on the 
harmonic oscillator during time r is given by Q 



P{W) 



5\W ~ (E'^ - E")]P'' P° 



(3) 



where = exp(— /3i?")/Zo is the initial (thermal) occu- 
pation probability and P^ „ are the transition probabil- 
ities between initial and final states n and m, 



pr 

m.7i 



dxo / dx(t)'^*{x)U{x,xo;T)(t)"^{xo) 



(4) 



The evaluation of the transition probabilities ^ requires 
the computation of the propagator U{x,xo',t) of the 
time-dependent harmonic oscillator. To this end, we fol- 
low the approach developed by Husimi [llj . 

The dynamics generated by the quadratic Hamiltonian 
^ is Gaussian for any uj{t). By introducing the Gaussian 
wave function ansatz, 



ipix^t) 



(^a{t)x^ +2b{t)x+c{t)'j 



(5) 



the Schrodinger equation for the quantum oscillator can 
be reduced to a system of three coupled differential equa- 
tions for the parameters a, b and c, 

^2 
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The nonlinear equation ([6]) is of the Riccati type. Writing 
a = it can be mapped to the equation of motion of 
a classical time-dependent linear oscillator, 



+ u;^{t)X = 



(9) 



The solutions of Eqs. ©-([H]), and hence the wave func- 
tion ([5]) of the quantum oscillator, can be recurrently 
constructed from the solution of Eq. ^ . An explicit so- 
lution of Eq. ^ for the linear parametrization |12l | , 



(10) 



is given in Appendix A. We stress that the following re- 
sults are valid for arbitrary uj{t). 

The general form of the propagator can be determined 
from the wave function ([5|) by noting that (l){x, r) = 
J dxU {x, xo',t) 4>{xo, 0). It is explicitly given by 



U{x,xo;t) 



2TTihX 



exp 



2nx 



{Xx^ -2xxo+Yxl) 



where X{t) and Y{t) are solutions of Eq. 

X{0) = 1(0) = 1 , 
y(0) = 1 Y{0) = . 



(11) 
satisfying 



(12) 



A compact expression for the generating function of the 
transition probabilities „ can then be obtained by us- 
ing the generating functions of the eigenfunctions 4>n{x) 
of the quantum harmonic oscillator. 



By exploiting the Gaussian character of Eqs. (Ilip and 
(Uni), the generating function of Eq. ^ becomes 



P{u,v) 



V2 



^Q*(l-M2)(l-t,2) + (l + „2)(l_|_„2)_4j, 



(14) 



where we have introduced the quantity 



Q* [ujl{u;lX{rr+X{T)') + {u^fY{Tr+Y{rr)] 

ZLUqUJi 

(15) 

The parameter Q* will play an important role in the fol- 
lowing discussion. We will indeed see that it can be re- 
garded as a measure of the degree of adiabaticity of the 
process. An interesting feature of the generating function 
P(u, v) is further that its (u, w)~dependence remains the 
same for all possible transformations uj{t). Details about 
the specific parametrization of the angular frequency only 
enter through different numerical values of Q*. 

With the help of Eq. ((T4)) . the characteristic function 
of the work, defined as the Fourier transform of the prob- 
ability distribution ([3]), 



G{fi) = / dWe'^"^ P{W) 



n=0 



oit) 



ri7r2(l-ti2) 



exp 



mLj{t) {1+u'^)(x^+Xq)—4uxxo 
h 2(l-«2) 



(13) 



can eventually be written in closed form in terms of the 
energies, Eq = hujQ, Si = huJi (Ae = Si — Eq) and the 
inverse temperature (3. We have 



V2(l - e-P^°)e 



g2ifj,ei ) (1 — g-2(i/J-l-/3)eo ) -|_ -|- g2^iei ^^J^ _|_ g-2(i/^+/3)eo ) _ 4 gif^si g-{il^+l3)eo 

I 



(17) 



The above expression for G(/x) is exact and fully charac- 
terizes the work distribution of the time-dependent quan- 
tum harmonic oscillator (U), for arbitrary parametriza- 
tion of the angular frequency uj(t). As mentioned pre- 
viously, different realizations of ult) will merely lead to 
different values of the parameter Q*. One should also 
note that the mass m of the oscillator does not explicitly 
appear in (flT)) . 

Recovering the Jarzynski equality As shown in 
Ref. ^] , the Jarzynski equality ([1]) can be recovered from 



the characteristic function G(/i), without having to com- 
pute its inverse Fourier transform, by simply putting 
fi — if]. In this way, we obtain from Eq. ([T7|) . 



sinh 
sinh 



=.-/3AF 



(18) 



The latter agrees with the direct derivation of the free 
energy difference of the quantum harmonic oscillator (see 
for instance Ref. p^). 
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III. WORK PROBABILITY DISTRIBUTION 

The direct analytic evaluation of the nonequilibriuni 
work distribution P{W) by Fourier inverting Eq. (fT7)) 
does not seem to be feasible in the general case. In this 
section, we derive approximate expressions of the work 
probability density for adiabatic and nonadiabatic pro- 
cesses, in the zero and high-temperature limit. 

Adiabatic case We base our discussion of adiabatic- 
ity on the equivalent classical harmonic oscillator ^ 
since the characteristic function G(u) is fully determined 
through its solutions X{t) and Y{t) For an adiabatic 
transformation, the action of the oscillator, given by the 
ratio of the energy to the angular frequency, is a time- 
independent constant. In other words, for a quasistatic 
process, we have the two adiabatic invariants. 



u{t) 



- ■ (19) 



From the definition (|15p of the parameter Q* , we see that 
in this case we simply have Q* — 1. As a consequence 
P{u,v) = (1 — uv)~^ and „ = 6m.n- The latter is 
an expression of the quantum adiabatic theorem: For in- 
finitely slow transformations, there are no transitions be- 
tween different quantum states. The characteristic func- 
tion (fT7|) accordingly reduces to 



G(m) 



(1 



(20) 



Equation (j20p further simplifies in the limit of low and 
high temperature. In the classical limit, hf3 ^ 1, we find 



G(m) 



(3ujQ — iAujfi 



(21) 



The inverse Fourier transform of Eq. (|2ip then yields the 
adiabatic work distribution 



P{W) 



e(w) 



(22) 



This result is identical to the classical work probability 
distribution derived by Jarzynski p^ . 

In the opposite limit of low temperature, h(3 3> 1, we 
have 



G(m) 



(23) 



As a result, the zero-temperature adiabatic work distri- 
bution is given by a delta function, 



P{W) = s(w 



(24) 



Nonadiabatic case We now turn to the nonadia- 
batic case Q* > 1. Wc start from the general expression 
(|17p for the characteristic function and expand it in the 
classical limit, <C 1. We obtain 



G(/i) 



y^Afl^ + B^l + C 



(25) 



where we have defined the constants A — 2Q*ujou>i—u>q — 
Lof, B = 2i(/3wg — Q* Puot^i) and G — /5^Wg. In order to 
perform the inverse Fourier transform of Eq. (|25p . we use 
the following integral [l^ . 



+ 00 



dx 



(x2 + z'^)P 



(26) 

where T[x) denotes the Euler Gamma function and 
K^{x) is the Macdonald function, i.e the modified Bessel 
function of the third kind. Combining Eqs. (|25p and (I26p . 
we then arrive at the nonadiabatic work distribution, 




W 



FIG. 1: Work probability distribution in the classical limit in 
the adiabatic and nonadiabatic case, Eqs. (|22|l and (|27|l . for 
the parameters too = 1, tJi = 1.3, (5 = 0.1 and Q* — 5. 



P{W) 



2/32^2/7r 



2Q*iLjouji 



rW 



2Q*L0{)Uii 
I 



\W\ 



(27) 



In the low temperature limit, h(3 ^ 1, the characteris- 
tic function pT|) reads, 

G{fi) ^ ^ . (28) 

^Q* + 1 - (Q* - l)e2v^i 



The inverse Fourier transform of Eq. (|28|) can be approx- 
imated in the limit of small £i by. 



P{W) 



^/tt{Q* - l)el(M^- Ae/2) 



(29) 
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FIG. 2: Work probability distribution at zero temperature in 
the adiabatic and nonadiabatic case, Eqs. (|24p and (|29|) . for 
the parameters eo = 0.1, £i = 0.13 and Q* = 5. 

The zero-temperature nonadiabatic work distribution 
(p9| is valid when IF > ^. 

IV. DISCUSSION 

The parameter Q* defined in Eq. (fT5|) can be ascribed a 
simple physical meaning . Let us consider a transition 
from initial state n to final state m. Then the average 
quantum number in the final state is 

(m)„=^mP;.^„ = (n+i)Q*-i , (30) 

rri 

which follows from the expression of the generating func- 
tion of Pmn- Similarly, the mean-square quantum 
number after time r is given by 

= ((m- (m)„)2)„ = i(g*2-l)(n2+n + l) . (31) 

Equation ((3T|) shows that the parameter Q* directly con- 
trols the magnitude of the variance cr^„ „. In the adiabatic 
limit, Q* = 1, we readily get (rn)„ = n and cr'^ „ = 0. 
We therefore recover that for quasistatic processes, the 
system remains in its initial state, m — n. On the other 
hand, for fast nonadiabatic processes, the mean (m)„ 
and the dispersion cr^ „ increase with increasing values 
of Q*, indicating that the quantum oscillator ends in a 
final state m which is farther and farther away from the 
initial state n. The behavior of Q* as a function of the in- 
verse switching time l/r is shown in Fig. (3) for the linear 
parametrization (|10p with fixed values of the initial and 
final angular frequencies. We observe, as expected, that 
infinitely slow switching, l/r —f 0, corresponds to the 
adiabatic limit Q* — > 1, while, for finite-time transfor- 
mations, the parameter Q* monotonically increases with 
the inverse switching time. 

The adiabatic and nonadiabatic work distributions 
Eqs. (EH), IMl) and Eqs. ^ are plotted in Figs. (1) 
and (2). They illustrate the effect of temperature and 



nonadiabaticity on the probability density P{W) of the 
harmonic oscillator. The mean and the variance of the 
work probability distribution, as well as the mean dis- 
sipated work can be found in Appendix B. Let us be- 
gin with Fig. (2). For an adiabatic transformation at 
zero temperature, the oscillator remains in the ground 
state for all times. As a consequence, the total work 
and the free energy difference are equal to the energy 
difference between the initial and final ground states, 
W = AF = Ae/2. In this case, P{W) is simply a delta 
function. As we now increase the temperature, the ini- 
tial thermal distribution broadens. This then results in 
a widening of the work distribution to a one-sided expo- 
nential in the high-temperature limit as seen in Fig. (1). 
On the other hand, when the nonadiabaticity parameter 
Q* is augmented, we observe still in Fig. (1) that the work 
distribution acquires a negative tail. This is a remark- 
able fact since for a process with negative total work, the 
quantum harmonic oscillator is cooled when its angular 
frequency is varied and not heated. This situation can 
obviously not occur when the oscillator is already in its 
ground state, explaining the absence of a negative tail 
in the nonadiabatic work distribution in Fig. (2). Here 
the broadening of the distribution P{W) is related to the 
work which is dissipated when the oscillator is brought 
to an excited state. 




2 4 6 e 10 12 14 1/x 



FIG. 3: Parameter Q* , given by Eq. (|15p . as a function of the 
inverse switching time 1 /r. Q* is here calculated for the linear 
parametrization (|10p with parameters ljo = 1 and ui = 13. 



V. SUMMARY 

We have analytically calculated the characteristic func- 
tion of the total work of a quantum harmonic oscillator 
with a time-dependent angular frequency. We have ob- 
tained approximate formulas for the work distribution 
for the case of adiabatic and nonadiabatic transforma- 
tions, in the limit of low and high temperature. We 
have further checked the validity of the quantum Jarzyn- 
ski equality for arbitrary parametrizations of the angular 
frequency and discussed the effect of temperature and 
nonadiabaticity of the work probability density. 
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In this appendix, we provide the analytical solutions of 
Eq. ([9]) for the linear parametrization (fTO| of the angular 
frequency, obeying the boundary conditions (|12p . They 
can expressed with the help of the Airy-A and Airy-B 
functions p^ . 



Appendix A 



The equation of motion for the time-dependent har- 
monic oscillator ([9]) can in general not be solved exactly. 



X{t) 



All 



Ai 



{t — t)lo1 — ttol 



Bi 



UJnT3 



[{uJo - wi)(a;o +cui)]3 



Bi 



{t — t)lUq — tuj\ 



T3 [(wo - Wl)(wo + t^l)]'' 



(32) 



and 



Y{t) 



TT Ai' 



Ai 



Bi 



T UJn 



-3 [{UJQ - UJl){uJli + Wl)] 3 



(t — t)uJq — tuj 



T3[(a;o -a;i)(wo +wi)]3 , 



\Bi'{ - 



LUnT'3 



[{ujo - wi)(wo +t^i)]3 



(33) 



where Ai'{x) and Bi'{x) denote the first derivatives with 
respect to x. 

Appendix B 

In this appendix, we evaluate the mean and the vari- 
ance of the work probability distribution, as well as the 
mean dissipated work for the frequency-dependent quan- 
tum harmonic oscillator. We use the formulas, 

{W) = / dWWP{W) = -iG'{0) , (34) 



We consider the general expression PT|) for the charac- 
teristic function of the work G'(/i) as well as the various 
approximations (I21|l . ((23|) . (|25|) and ([28|) obtained in the 
high and low temperature limit and for adiabatic and 
nonadiabatic transformations. The free energy difference 
AF is given by Eq. (fT2| . In the quantum limit S> 1, 
the free energy difference simplifies to 

AF = i(£i-£o), (37) 
whereas in the classical limit <C 1, one has 



AF^-ln — . (38) 
{W^)= dWW^P{W)^-G"{0), (35) P UJo 



Note that for an adiabatic transformation Q* = 1 and 
(Wdis) = (W) - AF . (36) Aoj/ujQ < 1, so that in this case AF ~ Aw//3. 
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Binh [-g-eil 


4sinh2 [feo] 
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General 
Q* = 1 
Eq. (HIl) 


i(£i - £o) coth [f eo] 


4(^1 ^0) si„h^[feo] 





?i/3< 1 

QVl 
Eq. (HID 


^(Q*wi - wo) 


_ 2Q*UJoUJi + (2Q*2 _ 1)^2) 


P \ UJQ ^UJQ 


1 

Q* = 1 
Eq. ([21]) 









n/3> 1 
Q* 7^ 1 
Eq. (IMl) 


i((3*ei - en) 


i(0*2_ 1)^2 
2 V ^ / 1 




?i/3> 1 
Q* = 1 

Eq. dSSl) 
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